load motor.mat;
XtrainNorm = motorNormalize(Xtrain);

X = linspace(-1, 1);
X = X'
YMeanPrediction = ones(size(X), 7);
Phi = zeros(size(X), 20);
for i = 1:31
   Phi(:, i) = X;
end

Phi(:, 1) = 1;
for j = 1:30
    %Phi(:, j+1) = Phi(:, j+1).^j;
    for x = 1:size(Phi(:, 1))
        Phi(x, j+1) = exp(-(Phi(x, j+1)-mu(j))^2 / 2*sigma^2);
    end
end

W = zeros(31, 31);

numDataPoints = size(XtrainNorm,1);
for Mindex=1:31
    powerIndex = Mindex;
    designMatrix = ones(numDataPoints,powerIndex);
    for rowDesign=1:numDataPoints
        for colDesign=2:powerIndex
            designMatrix(rowDesign,colDesign)=exp(-(XtrainNorm(rowDesign,1)-mu(colDesign-1))^2 / 2*sigma^2);
        end
    end
    w = (designMatrix' * designMatrix)\designMatrix' * Ytrain;
    w(31,1) = 0;
    W(:, Mindex) = w;
    
    if Mindex-1 == 5
        YMeanPrediction(:, 1) = Phi*w;
    end
    
    if Mindex-1 == 10
        YMeanPrediction(:, 2) = Phi*w;
    end

    if Mindex-1 == 15
        YMeanPrediction(:, 3) = Phi*w;
    end
    
    if Mindex-1 == 20
        YMeanPrediction(:, 4) = Phi*w;
    end
    
    if Mindex-1 == 25
        YMeanPrediction(:, 5) = Phi*w;
    end
    
    if Mindex-1 == 30
        YMeanPrediction(:, 6) = Phi*w;        
    end
    
    if Mindex-1 == 7
        YMeanPrediction(:, 7) = Phi*w;        
    end
end

figure(251)
subplot(2,3,1)
plot(X, YMeanPrediction(:, 1))
subplot(2,3,2)
plot(X, YMeanPrediction(:, 2))
subplot(2,3,3)
plot(X, YMeanPrediction(:, 3))
subplot(2,3,4)
plot(X, YMeanPrediction(:, 4))
subplot(2,3,5)
plot(X, YMeanPrediction(:, 5))
subplot(2,3,6)
plot(X, YMeanPrediction(:, 6))

figure(253)
plot(X, YMeanPrediction(:, 7))
